<config lang="json">
{
  "name": "UC2 Holorecon",
  "type": "native-python",
  "version": "0.0.2",
  "api_version": "0.1.2",
  "description": "Backpropagate an image into the sample plane.",
  "tags": [],
  "ui": [
          "Distance to the sensor plane: {id:'z_dist', type: 'number', min: 0, placeholder:1e-3}"
        ],
  "inputs": null,
  "outputs": null,
  "flags": [],
  "icon": null,
  "env": null,
  "runnable": false,
  "requirements": ["numpy", "matplotlib", "pillow", "scikit-image"],
  "dependencies": []
}
</config>



<script lang="python">
import numpy as np
import base64
import asyncio
import os
import time
from skimage import io
from PIL import Image
from io import BytesIO

import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
plt.ioff()

class ImJoyPlugin():
    def setup(self):
        self.window = None

    async def run(self, my):
        api.showStatus('[UC2] Python worker to reconstruct a hologram')
        print(my.config.chart_type)
        # read the raw-hologram
        my_holo_file = 'https://raw.githubusercontent.com/bionanoimaging/UC2-GIT/master/WORKSHOP/INLINE-HOLOGRAMM/PYTHON/hologram_mouse.jpg'
        myholo = 1.*io.imread(my_holo_file)
        z_dist = 1e-3
        await self.holorecon(myholo, z_dist)

    async def process_hologram(self, base64_url):
      ''' Bridging javascript with Python '''
      await api.showMessage('decoding image...')
      # remove encoding header
      data = base64_url.split('base64,')[1]
      # decode base64 to bytes
      img_bytes = base64.b64decode(data)
      buffer = BytesIO(img_bytes)
      # open the image as PIL image
      img = Image.open(buffer)
      arr = np.array(img)
      # process the image
      z_dist = 1e-3
      await self.holorecon(arr, z_dist)

    async def holorecon(self, myholo, z_dist):
        ''' start plugin here '''
        await api.showMessage('processing hologram...')
        # define some experimental parameters

        my_background_file = 'test2bg.jpg'
        mychannel = 2 # select the color channel 0,1,2
        mysize = 1024 # region of interest around the center coordinate
        # define acquisition parameters
        pixelsize = 1.4e-6; #0.6e-6#.2e-6#3.000e-6;
        lambda0 = 440e-9; # in nanometer #530e-9;
        # start and stop z-focus measure
        stepsize = 0.001;

        # select only the blue channel (because we used the blue LED)
        myholo = myholo[:,:,mychannel ]
        # crop out the ROI around the center to speed up computation (RADIX 2 please!)
        holosize = myholo.shape
        myholo =  myholo[int(holosize[0]/2-mysize/2):int(holosize[0]/2+mysize/2), int(holosize[0]/2-mysize/2):int(holosize[0]/2+mysize/2)]
        if(0):
            # display intermediate result
            plt.imshow(myholo)
            plt.colorbar()
            plt.show()
        # normalize hologram
        myholo = myholo/np.max(myholo);
        # estimate hologram's amplitude
        myamplitude = np.sqrt(myholo);
        ## Backpropagate the Hologram
        # Fresnel Propagation - test out a range of z-values to find the one that
        # brings the bug into focus
        Ef = FresnelPropagator(myholo, pixelsize, lambda0, z_dist)
        # Plot curve and save as png
        fig1, ax = plt.subplots(num='RECONSTRUCTED HOLOGRAM')
        plt.imshow(abssqr(Ef),cmap='gray')
        name_plot = 'HOLOGRAM.png'
        plt.savefig(name_plot,dpi=300)
        with open(name_plot, 'rb') as f:
            data = f.read()
            result = base64.b64encode(data).decode('ascii')
            imgurl = 'data:image/png;base64,' + result
            data = {"src": imgurl}
            data_plot = {
                'name':'Plot charts: show png',
                'type':'imjoy/image',
                'w':12, 'h':15,
                'data':data}
        ## Check if window was defined
        if self.window is None:
            self.window = await api.createWindow(data_plot)
            print(f'Window created')
        else:
            print(f'Update window.')
            try:
                await self.window.run(data=data)
            except:
                self.window = await api.createWindow(data_plot)
                print(f'Could not print to old window. New window created.')
        await api.showMessage('done!')

api.export(ImJoyPlugin())

def abssqr(x):
    # this is what a detector sees (only intensities)
    return np.real(x*np.conj(x))
def FT(x):
    # this only defines the correct fwd fourier transform including proper shift of the frequencies
    return np.fft.fftshift(np.fft.fft2(x)) # Fourier transform and shift
def iFT(x):
    # this only defines the correct inverse fourier transform including proper shift of the frequencies
    return np.fft.ifft2(np.fft.ifftshift(x)) # inverse Fourier transform and shift
def FresnelPropagator(E0, ps, lambda0, z):
    # Parameters: E0 - initial complex field in x-y source plane
    #             ps - pixel size in microns
    #             lambda0 - wavelength in nm
    #             z - z-value (distance from sensor to object)
    #             background - optional background image to divide out from
    #
    print('Entering Fresnelpropagator')
    upsample_scale = 1;                 # Scale by which to upsample image
    n = upsample_scale * E0.shape[1] # Image width in pixels (same as height)
    grid_size = ps * n;                 # Grid size in x-direction
    # Inverse space
    fx = np.linspace(-(n-1)/2*(1/grid_size), (n-1)/2*(1/grid_size), n)
    fy = np.linspace(-(n-1)/2*(1/grid_size), (n-1)/2*(1/grid_size), n)
    Fx, Fy = np.meshgrid(fx, fy)
    # Fresnel kernel / point spread function h = H(kx, ky)
    # from Fourier Optics, chapter 4
    # H = sqrt(z*lambda0)*exp(1i*pi*lambda0*z*(Fx.^2+Fy.^2));
    # sphere=exp(i*k/2/zc*(xx.^2+yy.^2));
    print('Generate Propagator.')
    print(str(Fx))
    print(str(Fy))
    print(str(lambda0))
    print(str(z))
    H = np.exp(1j*(2 * np.pi / lambda0) * z) * np.exp(1j * np.pi * lambda0 * z * (Fx**2 + Fy**2))
    #H= cos(pi*lambda0*z*(Fx.^2+Fy.^2)+(2*pi*z)/lambda0)+1i.*sin(pi*lambda0*z*(Fx.^2+Fy.^2)+(2*pi*z)/lambda0);
    # Compute FFT centered about 0
    print('Perform Fouriertransform.')
    E0fft = FT((E0));     # Centered about 0 since fx and fy centered about 0
    # Multiply spectrum with fresnel phase-factor
    print('Propagate.')
    G = H * E0fft
    print('Perform inverse Fouriertransform.')
    Ef = iFT(G) # Output after deshifting Fourier transform
    return Ef
</script>
